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Abstract 

An approach to land surface temperature (LST) estimation that relies upon 
Bayesian inference has been tested against multiband infrared radiometric 
imagery from the Terra MODIS instrument. The algorithm employed re- 
quires minimal knowledge of surface emissivity, starting from a parsimo- 
niously chosen (hence, uninformative) range of prior band emissivity knowl- 
edge. Two estimation methods have been tested. The first is the iterative 
contraction mapping of joint expectation values for LST and surface emis- 
sivity described in a previous paper. In the second method, the Bayesian 
algorithm is reformulated as a Maximum A-Posteriori (MAP) search for 
the maximum joint a-posteriori probability for LST, given observed sensor 
aperture radiances and a-priori probabilities for LST and emissivity. 

Two MODIS data granules each for daytime and nighttime were used 
for the comparison. The granules were chosen to be largely cloud-free, with 
limited vertical relief in those portions of the granules for which the sensor 
zenith angle \ZA\ < 30°. Level IB radiances were used to obtain 500 LST 
estimates per granule for comparison with the Level 2 MODIS LST product. 

The Bayesian LST estimators accurately reproduce standard MODIS 
product LST values. In particular, the mean discrepancy for the MAP re- 
trievals is I (AT) I < 0.3 i^, and its standard deviation does not exceed IK. 
The ±68% confidence intervals for individual LST estimates associated with 
assumed uncertainty in surface emissivity are of order 0.8 i^. 

The Appendix presents a proof of convergence of the iterative contrac- 
tion mapping algorithm. The expectation values of surface temperature in 
multiple bands, and jointly in all bands, converge to a fixed point, within a 
stipulated convergence criterion. Provided the support [Tjnin,Tmax] for the 
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calculation of the expectation value brackets the maximum in the joint pos- 
terior probability for LST, the fixed point coverges to the MAP LST estimate 
in the limit Tmax - Tmin 0. 



1. Introduction 



Land surface temperature (LST) is a vitally important remotely observ- 
able tracer of mass and energy exchange across the interface between the 
atmosphere and the ground. LST derived from satellite observations is of 
interest in its own right for local climate studies and climate change moni- 
toring, and as a component of studies of land co ver, land cover change, surface 



moist ure, and precision farming, among others ( iDash et a/.l . l2002l : IWan et al. 



20041 ) ■ Unfortunately, few surfaces at the bottom of the atmosphere radiate 
as blackbodies, and the parameterization of their surface state by means of 
emissivities complicates the task of accurate temperature determination. 

This paper continues the development of Bayesian LST estimators that 
do not require accurate knowl edge of surface emissivity, given radiance in 
multiple bands ( IMorganl . l2005l ). The approach may be considered as com- 
plementary to the widespread use of regression-law based split-window algo- 
rithms. 

Section [2] reviews the elements of the Bayesian LST formalism. Two im- 
plementations of the Bayesian approach have been investi gated. The first 



is the iterative contraction mapping algorithm presented in iMorganl (j2005[ ). 
The second is a Maximum-A Posteriori approach to LST estimation. Sec- 
tion [3] presents Bayesian estimates of LST and band emissivity for MODIS 
imagery. Section H] discusses results, and a final Section [5] offers conclusions. 
Appendix A presents a proof of convergence for the iterative algorithm, and 
estimates of the effect on LST retrieval of uncertainty in aerosol loading. 



2. Elements of Bayesian Land Surface Temperature Estimation 

2.1. The Posterior Probability for T 



We begin by recalling some of the results obtained in IMorganl (120051 ) . 
which may be consulted for details. The probability of sensing radiance /j 
at the top of the atmosphere (TOA), given T, e^, and an estimate o"j of the 
variance of the sensor radiance noise, is obtained by a Maximum Entropy 
(MAXENT) argument. The MAXENT estimator assumes the existence of 
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a forward model for the TOA radiance in terms of T and e^. For the case 
of a sensor aperture radiance forward model If^ that is linear in band 
emissivity e^, the posterior probability P{Ii \ T,ei,ai) is Gaussian in the 
mismatch between /,■ and If^^ in each band i. We also have 



N 



P{{I,} \T,{a,}) = l[P{I,\T,e„a, 



(1) 



1=1 



as the posterior probability for all bands jointly. 

Given radiance Jj in N bands detected at the top of the atmosphere (TOA) 
that originates from a patch on the Earth's surface by a sensor with noise 
radiance (Tj, the posterior probability that the surface is at a temperature T 
is given by Bayes' theorem as 



P(r,e, I I,,K) = P{T,e, \ K) 









K) 



(2) 



Equation ([2]) is evaluated with aid of the prior probability for the surface to 
be at temperature T and emissivity given available knowledge 



dT 

P(T, ei I K) dT dei oc — dei 



(3) 



The result is a posterior probability density. If one is ignorant of, or un- 
concerned with, the value of emissivity, it is possible to marginalize on the 
nuisance variable e, : 



P{T\h,K) 



deiP(T,e, I Ii,K) 



(4) 



and P{T \ li, Kj^ may b e expressed in c losed form using error functions. The 
explicit form is given in iMorganl (j2005[ ) . 



2.2. Land Surface Temperature Algorithms 

We desire to find an estimator for surface temperature using the poste- 
rior probability given by Bayes' Theo rem (|2l). Two LST algorithms have 
been developed. The first, recounted in iMorganI (120051 ) . iteratively computes 



^In lMorganl (j2005l ). the posterior probability density was written P{T \ Ii,K) 



dT 
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expectation values of LST and band emissivities from their posterior proba- 
bilities while contracting the range of a-priori limits on these parameters. 

The second calculational method is search for the maximum A-Posteriori 
probability (MAP) value of LST that maximises P{T \ {Ii},K). The moti- 
vation for investigating a MAP formulation for the Bayesian LST est i mator 



is simple. The original Bayesian LST algorithm presented in iMorganI (120051 ) 
obtains estimates of LST and band emissivity by iterative refinement of ex- 
pectation values for these quantities. Even though the posterior probability 
for surface T is a closed-form expression, this procedure is computationally 
intensive. Moreover, extensions to the estimator (such as, for example, a 
more careful treatmen t of th e forward model, or incorporation of calibration 



error effects (IMorganI . 120061 )) may deprive us of the comfort of the closed- 
form solution. The MAP criterion for an estimator promises to dramatically 
reduce the number of CPU cycles expended per LST estimate, and may yet 
preserve the ability to rely upon closed-form solutions. MAP LST estimates 
are obtained simply and quickly by a Go lden Section search]^ 



It was conjectured in IMorganI (120051 ) that the iterative algorithm con 



verges in general. Appendix A presents a proof that the iterative algorithm 
converges, and, in the limit, converges to the MAP estimate of LST. Granted 
the robustness of the Golden Section search, the existence of a solution for 
the MAP estimator may not appear to be a concern. However, the MAP 
algorithm in general includes iterative refinement of the a-priori limits on 
surface emissivity. The existence of a unique MAP estimate for LST in that 
case follows from a straighforward adaptation of the method of proof given 
in the Appendix. 

2.3. Confidence Intervals 

A desirable feature of an LST estimator is some means of assessing one's 
confidence in the accuracy of individual estimates. In the present case LST 
estimates are obtained by marginalizing over imperfectly known surface emis- 
sivities. It is naturally of great interest to explore means of judging the effect 
of an uninformative range for the a-priori probability for ej. Direct numeri- 
cal integration of the posterior probability P{T \ {Ii},K) with respect to T 
was used to calculate emissivity prior confidence intervals for LST retrievals 



^It is not difficult to extend the MAP approach to higher dimensions by simultaneously 
optimizing with Powell's method on (for example) band emissivity prior limits and LST. 
Sample tests of such extensions show promise, but remain to be investigated systematically. 
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presented in this study. Suppose that LST uncertainties arising from an un- 
informative emissivity prior are independent of those originating from other 
sources of error. Then the emissiviy prior confidence level may be combined 
with an estimate of the net effect of all other sources of error by adding in 
quadrature to obtain an overall LST uncertainty estimate. 



3. Application to MODIS imagery 

Both the MAP and iterative versions of the Bayesian algorithm have been 
applied to MODIS imagery. This effort took advantage of the availability of 
Level lb calibrated radiances and Level 2 land surface temperatures (used in 
lieu of ground truth) for Terra MODIS data granu les through the Goddar d 



Space Flight Center (GSFC) EOS/DAAC website fISavtchenko et aa 120041 ). 
Level 2 atmospheric profiles were obtained from the same site. The Level 

2 surface temperature product is used in preference to the standard Level 

3 product because it is in the same 1 km swath format as the Level 1 b 
radiance data. The GSFC-supplied HDFLook display package was used to 
manipulate the data granules and associated standard products. 

Four MODIS data granules were selected for comparison of Bayesian LST 
estimates with the M0D11_L2 LST product. To the extent possible, granules 
were chosen to have large areas with gentle topographic relief and minimal 
cloud cover effects. A list of the granules used in this study appears in Table 
1. Two granules are daytime, and two nighttime. The nighttime granule 
for 2006 day-of-year (DOY) 258 includes the alluvial plains of Iraq and both 
the coastal plains and central plateau of Saudi Arabia. 2006 DOY 268, also 
nighttime, includes the Great Plains portion of Central Canada. 2006 DOY 
350 is a daytime dataset covering Northeast Africa, including large plains in 
Sudan. 2006 DOY 347, also daytime, covers much of the plains of Central 
Asia. The M0D11_L2 product reported LST values for approximately 75% 
or more of the pixels in all four granules. 

Cloud-free pixels for which Level 2 MODIS LST values were reported 
were selected for analysis. A number of cuts were made in selecting data 
from the MODIS products for surface temperature estimation. The use of 
Level lb radiance products combined with Level 2 atmospheric profiles to 
generate LST for comparison with a Level 2 MODIS temperature product 
introduces sources of error originating in, among others, imperfect georef- 
erencing. Thus one may expect errors in surface altitude to affect retrieval 
accuracy. These effects may be exacerbated by use of bilinear interpolation 
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of 4 km horizontally reported profiles to the 1 km sampling of the spectral 
radiances. In order to minimize this source of potential error, LST retrieval 
was limited to sensor zenith angles \ZA\ < 30°. Cloudy pixels were masked 
by rejecting pixels for which the LST quality control mask bits 0,2,3,4, or 5 
were nonzero. LST retrievals were attempted only for pixels for which the 
dew point temperature at the base of the water vapor profile exceeded the 
freezing point of water. 

Cuts were also made on the basis of nonphysical or otherwise pathological 
fields in the granule data: Pixels for which the MODIS LST, the MODIS LST 
error, the surface temperature field in the profile granule (otherwise not used 
in the analysis), the surface air temperature at the lowest level of the profile, 
or the surface pressure were less than zero, were eliminated from processing. 
Finally, pixels for which surface pressure exceeded 1000 millibars, but for 
which the surface altitude was less than zero, were omitted from the analysis. 

Neither ground truth nor surface meterological range (VSBY) is avail- 
able as a standard product with MODIS granules. However, the presence 
of LWIR ban ds in the split-window algor i thm used to produce the MODIS 
LST product (IWan and Dozierl . ll996l : IWanl . ll999l ) suggests that LST retrievals 
should not be too sensitive to the exact value assumed for VSBY in the for- 
ward model. Accordingly, this study replicates MODIS LST values assuming 
a nominal, but fairly large, value of 23 km for the surface VSBY parameter 
for clear pixels. Appendix B presents a calculation which supports the con- 
tention that, for purposes of the present study, this approximation should 
not introduce large errors for clear air pixels. 

A random selection of 500 pixels was made from each granule, subject 
to data cuts described above. LST retrievals from both single-pass MAP 
and two-pass iterative versions of the algorithm were obtained for all gran- 
ules. In addition, ±68% confidence interval^ were computed for the MAP 
retrievals by numerical integration of the posterior probability distribution. 
The handset used in the retrievals appears in Table 2. T he MODerate resolu- 
tion atmospheric TRANsmission (M0DTRAN4) code (IKniezys et a/.l . Il996 : 



Berk et al. 



m 



1989) was used for the forward model calculations, as described 



Morgan! (120051 ). The forward models incorporated MOD07_L2 air tem- 



•^±68% confidence levels for lack of emissivity knowledge were chosen to approximate 
the equivalent of 2 cr intervals, given that P{T \ {li}, K) generally displays some skewness. 
This should be borne in mind when when comparing with the effect of other error sources 
described by a standard deviation. 
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perature and water vapor profiles and assumed a CO2 mixing ratio of 382.5 
ppm. 

Initially, the prior probability in (|3]) was assumed to hold over the range 

200/s: < T < 500/s: 
0.85 < e, < 0.999 

for T and e in each band. However, it proved advantageous to restrict the 
range of the emissivity prior for MODIS bands 31 and 32, 

0.95 < egi < 0.999 
0.95 < £32 < 0.999 

so as to take advantage of the observation in (Wan and Dozier 1996, Wan 
1999) that most natural land cover types have an emissivity in excess of 0.97 
in these bands. 



In the simulated LST retrievals presented in iMorganl ( 120051 ). the equiv- 
alent noise radiance cxj in each band was parameterized by a signal-to-noise 
ratio (SNR) value obtained by approximately inverting the noise-equivalent 
temperature change (NEAT) for the specified nominal MODIS performance. 
Experimentation with actual MODIS data revealed that the SNR values so 
chosen greatly underestimated the random variance in radiance measured by 
the actual Terra spacecraft against real terrestrial surfaces, through a real, 
and sometimes dirty, atmosphere. In consequence, the SNR parameter was 
reduced from values characteristic of ideal performance, in which the only 
source of radiance error in sensor measurements is accounted for by NEAT, 
by a full order of magnitude. In fact, the label "SNR" is probably somewhat 
misleading, since the random variability it describes is likely the concate- 
nation of numerous confounding factors, not all of which need accurately 
described as "noise". It is retained because of its use in earlier studies. 

The SNR values used in all retrievals reported in this study appear in 
Table 2. These values were selected on the basis of trial- and-error with 
pixels from DOY 347 that were drawn independently from the selection used 
for comparison with the MODIS LST product. LST estimates are fairly 
insensitive to the exact values chosen for band SNR's, so long as these are 
fairly large. If the SNR values are too large, the noise variance estimate 
becomes so tight that the retrieval becomes pathologically overconstrained. 
In the most commonly observed pathology, one or more of the individual band 
posterior probability fails to overlap the others; i. e. the support of the joint 
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posterior probability becomes the null set. If SNR is too small, < 0(10) for 
MODIS, LST retrievals become badly biased towards underestimates, and 
confidence intervals become worrisomely large. 

Tables 3 and 4 present the mean, standard deviation and of the mis- 
match for each granule for MAP and iterative retrievals, respectively. In 
Table 3, summary statistics for LST are shown for 0.85 < €31^32 < 0.999 as 
well as for 0.95 < €31 32 < 0.999. A small number of pixels required relaxation 
of the SNR parameter or emissivity priors in orde r to obtain an it erative LST 



estimate, in a pattern similar to that reported in iMorganl ( l2005l ). 

Figures 1-8 display plots of the mismatch between MAP and MODIS 
LST values, with ±68% confidence intervals, and of mismatch histograms, 
for 500 MODIS pixels from each of the four data granules. Overplotted on 
each mismatch histogram is the histogram for the equivalent Gaussian error 
distribution 



Note that the confidence intervals in the plots (whose mean and standard 
deviation for each granule is given in the LST mismatch figure captions) 
are for the uncertainty in LST resulting from the assumed uncertainty prior 
for surface emissivity only. Thus, they do not (directly) reflect sensor noise 
effects or other sources of uncontrollable variance]^ 

4. Discussion 

Agreement between the Level 2 MODIS LST product and MAP LST es- 
timates obtained in this study is better than 1 i^' in both mean and standard 
deviation of the mismatch for all granules. On average, the mismatch with 
MODIS LST in Table 3 increases only slightly for the looser emissivity prior. 
In addition to the tendency for the Bayes LST retrievals to display a slight 
negative bias with respect t o the MODIS product (consistent with the ob- 
servation made in [Morgan the mismatch plots show evidence of a 



positive trend in the discrepancy with respect to MODIS LST. The trend is 



^These other sources of variance do enter the confidence interval through the MAXENT 
estimator for the mismatch between sensor and forward model radiance. Thus, the width 
of the confidence interval for emissivity uncertainty depends upon what value for the SNR 
parameter is used in the retrieval. 
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readily apparent to the eye for the daytime retrievals, but can also be dis- 
cerned in Figure 3 for (nighttime) DOY 268. MAP and iterative solutions 
reproduce MODIS LST comparably well. Only the DOY 347 retrieval shows 
a notably different mean mismatch compared to MODIS between iterative 
and MAP approaches. For both algorithms, per degree of freedom [DOF) 
is close to unity. 

The best estimate of MODIS LST accuracy lies in the range ±0.6 A' to 
±1.0 K flWan et all \20Q4 . ICuenther 6^011 . l2002h . The standard deviation of 
the MAP temperature mismatch in Table 3 also lies in this range. A t-test 
shows the difference in means between MODIS and MAP LST values to be 
statistically significant for both nighttime granules and marginally (at the 
5% level) for DOY 347. However, the generally small level of the discrep- 
ancy indicates that, given better control of error sources, the Bayesian MAP 
algorithm is capable of replicating the accuracy of the MODIS M0D11_L2 
product. 

A small number of retrievals for both algorithms have large discrepancies 
with respect to M0D11_L2 LST. The cosmetic effect of these retrievals in 
the mismatch plots may give an exaggerated impression of their effect on the 
accord between MODIS and Bayesian LST. Judging from the mismatch his- 
tograms, they appear to be outliers resulting from pixels whose radiometric 
accuracy has been compromised by one or more of the error sources described 
earlier, despite imposition of the data cuts. In the case of the granule for 
DOY 2006/350, two pixels out of 500 have a discrepancy in excess of AK. 
These discrepant values appear significantly out-of-family with respect to the 
largely Gaussian appearance and favorable statistic for the retrievals as a 
whole. In Table 3, the summary statistics for DOY 350 are replicated with 
the two worst outlier pixels removed in the rows marked with an asterisk, 
with minimal change to the results. 

Error sources arising from georeferencing errors in Level lb radiance prod- 
ucts have already been discussed. An additional possible origin for the dis- 
crepancies with MODIS LST is errors in the forward model. Among other 
causes, these could result from the fixed choice of VSBY = 23 km {vide. Ap- 
pendix B, however), or occur as a consequence of using erroneous values of 
air temperature or water vapor from the MOD07X2 profiles in MODTRAN 
(Wan et al. 2004). Another posible source is imperfect cloud masking in the 
Level 2 MODIS products that could corrupt pixel radiances. This last source 
of error would afflict MODIS LST as well, and would have to affect MODIS 
and Bayesian LST differentially in order to contribute to discrepancies. Fi- 
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nally, the assumed band emissivity prior range need not always bracket the 
actual surface emissivity for a pixel. 

A noteworthy feature of the LST retrievals reported for this study is the 
comparatively small uncertainty for individual LST estimates originating in 
the assumed lack of surface emissivity knowledge. The mean value of the 
±68% interval for individual granules is greatest for DOY 350 (daytime) at 
0.828 ± 0.055 K and least for DOY 268 (nighttime) at 0.638 ± O.OUK. As 
noted earlier, the emissivity uncertainty confidence interval depends upon the 
assumed level of random variance due to noise and other uncontoUable factors 
capable of influencing LST retrieval. However, insensitivity of Bayesian LST 
to emissivity seems likely to hold generally. In particular, confidence intervals 
for €31^32 > 0.85 differ little from those for 631 32 > 0.95 

It was remarked earlier that the Bayesian approach to LST estimation 
was complementary to the more common split-window algorithmic family. 
The MODIS generalized split-window algorithm that furnished LST values 
against which performance of the Bayes approach has been tested relies upon 
regression against a large training dataset, and requires assignment of pixels 
to emissivity classes. The Bayesian approach, while statistical in its treat- 
ment of emissivity uncertainty, relies upon radiative transfer theory, rather 
than regression of brightness temperatures against true LST, to invert radi- 
ance to surface temperature. It is thus a method of remotely determining 
LST that (so to speak) minimizes hostages to the fortunes of emissivity 
knowledge. The accurate reproduction in the mean of M0D11_L2 LST in 
this study suggests the possibility of using the Bayesian approach for virtual 
"ground truth" to train local or regional split- window LST algorithms. 

Although the golden section search for the MAP version runs faster than 
the iterative solution in protoype IDL (Interactive Data Language) code, the 
difference in execution time is not great, and is in any event dominated by the 
execution time for the forward model calculations in MODTRAN. However, 
the MAP version of the algorithm would seem preferable both on the basis 
of simplicity and because it readily admits extension to higher-dimensional 
searches for optimal joint estimates of LST and band emissivity. 

5. Conclusion 

The performance of the Bayesian land surface temperature retrieval al- 
gorithms has been tested in a small-scale study by comparison of Bayesian 
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LST estimates against the MODIS M0D11_L2 LST product. Both the orig- 
inal iterative contraction mapping search and a new Maximum A-Posteriori 
algorithm reproduce MODIS LST values with good accuracy. The mean 
mismatch with respect to the M0D11_L2 LST product for Maximum A- 
Posteriori estimates does not exceed 0.3 K, and its standard deviation is less 
than 1 K, for any granule used in this study. For the iterative algorithm, 
the corresponding mean mismatch does not exceed 0.5 i^', nor its standard 
deviation, lAK. 

The agreement of the MAP estimates derived from MODIS Level IB 
radiance with the M0D11_L2 LST is comparable to the claimed accuracy 
for MODIS LST. Emissivity prior confidence intervals for individual LST 
estimates shown in Section H] show that LST retrieval accuracy with the 
Bayesian approach is insensitive to modest uncertainties in surface emissivity. 
We conclude the MAP LST algorithm is a viable candidate for LST retrieval 
in circumstances in which band emissivities are poorly known. 

It has been shown {vide. Appendix A) that the iterative r efinemen t of ex- 



pectation values underlying the TES algorithm presented in iMorganI ([20051) 
converges, in general. In the limit of vanishing support for the LST expecta- 
tion values, the iterative approach is equivalent to the Maximum A-Posteriori 
estimator for LST. 
Acknowledgement: 
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Appendix A. Convergence of iterative LST retrieval 



The algorithm used in IMorganI (120051 ) iteratively refines the calculation of 



the expectation values (IA.ip - (IA.3P below of surface temperature and emissiv- 
ity in multiple bands by systematic contraction of the limits of integration. 
It has been found that this procedure, which is described in detail in (Mor- 
gan 2005) , converges rapidly and reliably in practice. However, to date there 
has been no proof that the algorithm actually converges in general. This 
Appendix supplies such a proof. 

For reference, the expectation values for T given radiance in bands i, 
obtained from the joint posterior probability for observing radiances Ii,i = 
1 , N are 



(T,) = , (A.l) 

ff--P{T\I„a.) 

J. m in 
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assuming T is known to lie between a minimum and a maximum, while a 
joint estimator for T given radiances in all N bands is 



J^--P{T\{h},{a,}) 



An estimator for the emissivity in band i is given by 



\^'V rtmax r>//rr\ \ T \J K-^-'-'J 

J , P{{Ti},e Ii,ai)de 



Appendix A. 1. Assumptions 

Both surface temperature and band emissivities are assumed to lie within 
limits imposed a-priori. The surface temperature is limited to the range 

and emissivity to the range 

(^min ^ ^ ^max (-^■^) 

where the minimum and maximum values may be band-dependent. 

It is assumed that the maximum of P{T \ {Ii\,K) is always included 
in the integration. In any workable procedure for iteratively refining surface 
temperature or emissivity by contracting the support [T^j^, Tmax] of the inte- 
grations (]A.mA.3l) which define (T) and (ej), that support must bracket the 
maximum of P{T \ {Ii],K) as it contracts. Otherwise, the most probable 
values of (T) and (cj) would be excluded from the search, in violation of the 
MAXENT assumption that expected values should maximize the likelihood. 

Appendix A. 2. Proof of Convergence 

Proposition: The iterative refinement of (]A.ip - flA.2p converges. 



For clarity we limit the discussion to (lA.ip and (]A.2p . The extension to 
accomodate (1A.3P is clear. A total of iV -|- 1 expectation values for surface 
temperature (Tj) (the index i now including the joint estimate) is obtained 
at each stage of iteration. The values (Tj) comprise a vector in R^"*"^. These 
are confined within the subset of R^+^ delimited by flA.4p that provides 



nonvanishing support for the integrations in the definitions (lA.ip and f lA.2p . 
Call the vector of expectation values (T). 
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At the n-th stage of iteration, (T) is calculated over the range T^^^] 
and the mismatch amongst the band LST expectation values (Tj) is scored 
by 

Ai{T)) ^ maxm - {T,)\ (A.6) 

taken over all pairs of estimates. Let the set X = {x} be the set of all 
m-tuplets of the form 

— I max min \ ( A 7\ 

V m)-{T,)\,^^,J<N + l} ) ^^-'^ 
consistent with the overall limits (lA.4p . Here 

+ (A.8) 

where the quantity in parentheses is the binomial coefficient. The set X is 
compact, so all components of any sequence G X are bounded. (They are all 
nonnegative, as well.) Given points x^,x^ e X define a metric function 

d(x\x2)= (A.9) 



The p air (X, d) can be shown to be a complete metric space (lAliprantis and Burkinshawl . 



19811 . p. 35). 



Define the function : X — t- R by 

C(x) = eYxi>Q. (A.IO) 

i=l,Af 

where 9 > 1. 

Consider a contractive sequence in X with n = 1, 2, ■ ■ ■ for which 

> Vn (A.ll) 

Define the mapping from one iterate of x to another in X by 

x"+i = M(x"). (A.12) 

The function : X — )■ R is lower semicontinuo us, and M : X — )■ X is 
continuous, on X. By Caristi's fixed-point theorem (jOkl . 120071 ) the mapping 
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M has a fixed point M(x) = x. The iterative refinement of (T) therefore 
converges in X. Moreover, if 



hm = T°° 

mm mm 

n— ^oo 



then 



hm T^,. = T~,, (A.13) 
< (T,) < TZx (A.14) 



by the definition of (T). □ 

Appendix A. 3. Iterative algorithm and maximum a-Posteriori probability search 

The fixed point of the iterative mapping drives A — )• A°°, returning an 
N + 1-tuple of LST expectation values, none of which differ from any of the 
others by an amount greater than A°°. It may appear that nothing in the 
result just obtained says any such estimate necessarily gives an accurate value 
for LST. In practice, however, it happens that a small number of iterations 
suffices to yield accurate LST estimates. 

The manner in which T^ax — T^m contracts remains at our disposal. We 
are at liberty to force it to tend to zero. We now show that T^ax — Tmm — ^ 
is equivalent to a contractive mapping chosen with 

> x1+\ Vn (A. 15) 



rather than flA.lip . and that as — > oo (Tj) converges to the Maximum 



A-Posteriori estimate Tmap- 

The support for the computation of (T) tends to zero as a; — > 0. There 
is thus a nested sequence of intervals [T"1\T" ^J] C [T".^,r"„^]. Cantor's 



theorem on the intersection of a nested sequence (lAlipr ant is and Burkinshawl . 



19811 . p. 32) implies the limit of this process will be a singleton in X. It was 
assumed at the outset that the support of integrations [T^j„,T^„^] in ( lA.lt - 
IA.2P brackets Tmap', the singleton thus corresponds to T = Tmap, which is 
also the value of (Tj). We may see this last by the following argument: In 

/J™^TP(r I {I^},{(ri}) 
J^--P{T\{h},{a,}) 
as Tynin,Tynax Tmap, wc may write 

/oo 
dTf{T)T (A.17) 
-oo 
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with 



f{T) 



P{T\{h},{^.}) 

else 



T < T < T 

; mm _ _ -"^ max 



(A.18) 



As Tjnax — Train — )■ 0, it caii be shown ( iLighthilll . Il960l ) that 



f{T)^6{T-TMAp) (A.19) 

and thus 

(Tj) ^ Tmap. (A.20) 

As X — )■ 0, A — )■ and the other expectation values (Tj) — > Tmap, a 
unique converged surface temperature estimate. Conversely, if T^^^—T^^^ — >■ 
0, (T) — )• Tmap by applying the preceding argument to each expectation 
value separately, and x — )■ 0. Then a subsequence x"' exists that respects 
(lA.lSp . We thus have the 

Corollary: In the limit T^^x ~ ^mm 0; iterative refinement of unbiased 
expectation values, i.e., those for which supp{T) brackets the maximum of 
the P{T I {Ii},K), is equivalent to the Maximum A-Posteriori (MAP) ap- 
proach, in which one searches directly for the value of T which maximizes 
the joint posterior probability ([1]). □ 

This result is likewise true of the straightforward extension of the proof 
given here to the case of a contractive iterative search for a fixed point for 
(T,e.). 



Appendix B. Insensitivity of LST estimates to meteorological range 



The contribution of aerosol loading of the atmosphere to the forward 
model radiance calculations was parameterized by specifying the surface me- 
teorological range, or visibility (VSBY). This is the variable used by MOD- 
TRAN to scale aerosol extinction. The meteorological range was assumed 
to be 23 km for all retrievals presented in the main body. Because LST re- 
trievals were only attempted for cloud- free pixels in what were judged to 
be reasonably clear air conditions, it seems reasonable to suppose that this 
assumption should not introduce significant error so long as few pixels have 
VSBY < 23 km. We now examine the likel y validity of th is assu mption. 

We begin by adapting the analysis of iMorganI (120051 . 120061 ) to obtain 
the prior probabilty for surface meteorogical range VSBY and use it to 
construct posterior probabilities marginalized over a range of VSBY. The 
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starting point is the apparently trivial observation that the meteorological 
range is a length: Traditionally, and loosely, VSBY is the distance at which 
an observer looking horizontally at a height of two meters just fails to discern 
the contrast presented by an object against the horizon. In the remote sensing 
of objects by electromagnetic radiation, Poincare invariance requires that any 
two observers must be able to relate their description of events by a Lorentz 
transformation. In particular, the relation between VSBY = V as observed 
at time t in primed and unprimed coordinate systems moving with relative 
velocity v = (3 c with respect to each other along a stipulated line of sight 
(which might as well be taken to be parallel to the ground plane) is given by 

V^^{V'-^t') (B.l) 

t = 7(i' - ^ V') (B.2) 

where 

7 = , ^ . (B.3) 

The two observers will assign prior probabilities for the occurence of par- 
ticular values of the meteorological range 

f{V)dV (B.4) 

in the unprimed system and 

g{V')dV' (B.5) 

in the primed one. In order for observers to agree as to the form of the 
estimator, the priors must be related by 

g{V')dV' = J-^f{V)dV (B.6) 

where 

dV 

is the Jacobian for the transformation between descriptions in the parameter 
space. Primed and unprimed observers, if equally cogent, must also agree as 
to the functional form of the prior probability for V ^ thus 

m = m (B.8) 



16 



leading to a form of Schroder's equation for /, 



f{V)=jf{jV), (B.9) 

with solution 

nv) = ^ (B.io) 

for the prior probability. We conclude that the VSBY parameter, with 
dimensions of a length, is a sca le para n aeter, for which the prior probability 



takes the Jeffreys form fIB.lOl) flJayned . Il968l ) 



The prior probability f lB.lOp may be inserted into an estimator for LST 
which marginalizes on both surface emissivity and VSBY in each band i, 



VSBYy 

VSBY< 



lax -J T/" 

P{T\h,a,)= I I P{T\h,e,V,a,)de—, (B.ll) 



from which a joint posterior probability may be calculated as a function of 
LST. This quantity may be used to compute MAP estimates of LST and con- 
fidence intervals. The estimator will now be used to examine the sensitivity 
of LST estimates when VSBY was allowed to vary over a range of values by 
means of an illustrative sample calculation. 

LST MAP retrievals and confidence intervals as a function of VSBY^ for 
one of the pixels used in the DOY 350 granule comparison with MODIS LST 
appear in Figure The ±68% confidence intervals shown account for the 
effects of the VSBY and emissivity priors only, and include no estimate of 
uncertainty originating from noise effects. The upper limit on the visibility 
prior was fixed at VSBYy = 50 km. The surface T prior for these retrievals 
was constrained to the range 306 K <T< 317 K. 

The values of LST and associated confidence intervals for this example 
cannot be directly compared with the MAP value 311.79 K marginalized on 
surface emissivity alone. However, note that the range of LST variation in 
Figure 9 is limited, even within the constrained LST prior limits. Moreover, 
the confidence intervals tend to a fixed range (consistent with the MODIS 
LST value) as VSBY increases. These observations provide a degree of 
empirical support for the assumption that the LST retrievals presented in 
the main text should not be too sensitive to the exact value of VSBY used 



calculation for a single pixel is given because the use of nested Romberg quadratures, 
in which the innermost loop invokes MODTRAN, is computationally intensive. 
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to parameterize aerosol loading in the forward model, so long as the actual 
VSBY is great enough. 
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Table 1. MODIS granule datasets 



MOD021KM.A2006258. 1900. 005. 2006260185940. hdf 
MOD021KM.A2006268.0455.005.2006270054543.hdf 
MOD021KM.A2006347.0615.005.2006348125833.hdf 
MOD021KM.A2006350.0830.005.2006352011824.hdf 



MODll_L2.A2006258.1900.004.2006259185248.hdf 
MODll_L2.A2006268.0455.004.2006269162220.hdf 
MODll_L2.A2006347.0615.004.2006349135340.hdf 
MODll_L2.A2006350.0830.004.2006351140447.hdf 



MOD07_L2.A2006258.1900.005.2006260192510.hdf 

MOD07_L2.A2006268.0455.005. 2006270062431. hdf 
MOD07_L2.A2006347.0615.005. 2006348132045. hdf 
MOD07_L2.A2006350.0830.005.2006352015143.hdf 



Table 2. MODIS band definitions, SNR parameters 
MODIS band wavelength hmits SNR 



20 3.660-3.840 /xm 25 

22 3.929-3.989 nm 25 

23 4.020-4.080 /xm 25 
29 8.400-8.700 /xm 25 

31 10.870-11.280 nm 50 

32 11.770-12.270 nm 50 
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Table 3. Maximum A-Posteriori surface temperature results: mean mis- 
match, standard deviation, maximum mismatch, and P^r degree of free- 
dom (DOF). Results are given for two choices for the lower limit to the 
emissivity prior in bands 31 and 32. The rows labeled with an asterisk dif- 
fers from those immediately above it by omission of the pixels with the worst 
two mismatches compared to MODIS LST 

Granule {LST mismatch) \Max mismatch] ^jDOF 



0.95 < €31,32 < 0.999: 






258.1900 


0.238 ± 0.723 X 


3.121 X 


1.106 


268.0455 


-0.286 ±0.416 


3.039 


1.470 


347.0615 


-0.058 ±0.672 


2.790 


1.005 


350.0830 


-0.036 ± 1.043 


6.974 


0.999 


350.0830* 


-0.059 ± 0.975 


3.462 


1.002 


0.85 < €31,32 < 0.999: 






258.1900 


0.454 ± 0.738 K 


3.082 K 


1.376 


268.0455 


-0.209 ± 0.429 


2.980 


1.236 


347.0615 


0.100 ± 0.731 


3.249 


1.017 


350.0830 


0.216 ± 1.109 


8.866 


1.036 


350.0830* 


0.189 ±1.019 


3.737 


1.032 



Table 4. Iterative surface temperature results: mean mismatch, standard 
deviation, maximum mismatch, and jDOF 

Granule {LST mismatch) \M ax mismatch] ^jDOF 

258.1900 0.263 ±0.690 2.856 X 1.143 

268.0455 -0.295 ±0.425 3.220 1.480 

347.0615 -0.670 ±0.966 4.081 1.480 

350.0830 0.038 ±1.355 7.789 0.999 
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Figure .1: Mismatch between MAP and MODIS LST for 500 nighttime pixels from DOY 
2006/258 granule. Mean confidence interval (AT) = 0.730±0.028 if . Confidence intervals 
reflect effects of emissivity uncertainty only. 
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Figure .2: Nighttime MAP vs. MODIS LST mismatch histogram for DOY 2006/258. Solid 
line: MAP-MODIS mismatch; dashed line: Equivalent Gaussian histogram corresponding 
to (ST) = 0.238K, a = 0.72iK. 22 
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Figure .3: Mismatch between MAP and MODIS LST for 500 nighttime pixels from DOY 
2006/268 granule. Mean confidence interval (AT) = 0.638 ± 0.014 if 
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Figure .4: Nighttime MAP vs. MODIS LST mismatch histogram for DOY 2006/268. 
Equivalent Gaussian histogram: {5T) = -0.286ii', a = 0A16K. 
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Figure .5: Mismatch between MAP and MODIS LST for 500 daytime pixels from DOY 
2006/347 granule. Mean confidence interval (AT) = 0.796 ± 0.054 if 
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Figure .6: Daytime MAP vs. MODIS LST mismatch histogram for DOY 2006/347. Equiv- 
alent Gaussian histogram: {ST) = -OMSK, a = 0.672ii'. 
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Figure .7: Mismatch between MAP and MODIS LST for 500 daytime pixels from DOY 
2006/350 granule. Mean confidence interval (AT) = 0.828 ± 0.055 iiT 
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Figure .8: Daytime MAP vs. MODIS LST mismatch histogram for DOY 2006/350. Equiv- 
alent Gaussian histogram: {ST) = -OmGK, a = 1.043K. 
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Figure .9: MAP LST estimates, confidence intervals marginalized on emissivity and VSBY 
vs. VSBY< for DOY 2006/350. The horizontal lines show the MODIS LST value 311.66± 
0.6K for this pixel. 
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